#include "vol.hpp"
#include <cstdlib>
#include <sstream>
#include <queue>
#include <vector>
#include <tr1/unordered_map>
#include <lpt/lpt.hpp>
#include <mytl/comb.hpp>
#include <mytl/point.hpp>

typedef lpt<3,lpt3_64_bits> lpt3d;
typedef lpt_tree<lpt3d> lpt3d_tree;

using namespace std;
using namespace mytl;

namespace std { namespace tr1 {

  template <>
  struct hash<point<3,double> > {
    size_t operator()(const point<3,double>& p) const {
		  unsigned int x=(p[0]+1.0)/2.0*1024;
		  unsigned int y=(p[1]+1.0)/2.0*1024;
		  unsigned int z=(p[2]+1.0)/2.0*1024;
			morton3_10(x);
			morton3_10(y);
			morton3_10(z);
      x=x | (y << 1) | (z << 2);
			return static_cast<size_t>(x);
  	}
  };

} }

double error(lpt3d s, norm_vol * v, int m) {
	array<point<3,double>,4> vs;
	array<vec<3,double>,4> ns;
	s.simplex(&vs[0][0]);
	for(int i=0;i<4;++i) {
	  vs[i][0]=(vs[i][0]+1.0)/2;
	  vs[i][1]=(vs[i][1]+1.0)/2;
	  vs[i][2]=(vs[i][2]+1.0)/2;
	  ns[i]=v->grad(vs[i][0],vs[i][1],vs[i][2]);
	}
	vec<4,double> mw=_((double)m,0.0,0.0,0.0);
	double meant=0.0;
	int count=0;
	do {
	  point<4,double> w=(1.0/m)*mw;
		point<3,double> p=barycentric_combination(w,vs);
	 	vec<3,double> g=v->grad(p[0],p[1],p[2]);
		normalize(g);
	 	vec<3,double> ig=barycentric_combination(w,ns);
		normalize(ig);
	  double c=dot(ig,g);
		double t=sqrt(1.0/(c*c)-1.0);
		meant+=t;
		count++;
	} while(next_triangular(m,&mw[0],&mw[0]+4));
	meant/=count;
	meant*=1.0/6.0/(1<<s.simplex_level());
	return meant;
}

int main(int argc, char * argv[]) {
  vector<string> arg(argv,argv+argc);
	if(argc<2) {
	  cerr<<"usage: vol2mesh [options] input_file.raw [output_file]"<<endl;
	  cerr<<"options:"<<endl;
	  cerr<<"\t --dim-x nx"<<endl;
	  cerr<<"\t --dim-y ny"<<endl;
	  cerr<<"\t --dim-z nz"<<endl;
	  cerr<<"\t --dim n"<<endl;
	  cerr<<"\t --type uchar|char|ushort|short|uint|int|float|double"<<endl;
	  cerr<<"\t --max-level Ml"<<endl;
	  cerr<<"\t --min-level ml"<<endl;
	  cerr<<"\t --max-error Me"<<endl;
		exit(1);
	} 
	string in_file="";
	string out_file="";
	string type="uchar";
	int dim_x=64;
	int dim_y=64;  
	int dim_z=64;  
	int min_level=2;
	int max_level=6;
	double max_error=0.00005;
	for(int i=1;i<argc;++i) {
	  if(arg[i]=="--dim") {
		  stringstream(arg[++i])>>dim_x;
			dim_z=dim_y=dim_x;
		} else if(arg[i]=="--dim-x") {
		  stringstream(arg[++i])>>dim_x;
		} else if(arg[i]=="--dim-y") {
		  stringstream(arg[++i])>>dim_y;
		} else if(arg[i]=="--dim-z") {
		  stringstream(arg[++i])>>dim_z;
		} else if(arg[i]=="--type") {
		  type=arg[++i];
		} else if(arg[i]=="--min-level") {
		  stringstream(arg[++i])>>min_level;
		} else if(arg[i]=="--max-level") {
		  stringstream(arg[++i])>>max_level;
		} else if(arg[i]=="--max-error") {
		  stringstream(arg[++i])>>max_error;
		} else {
		  if(in_file=="") in_file=arg[i];
			else out_file=arg[i];
		}
	}
	if(in_file=="") {
	  cerr<<"Missing input file."<<endl;
		exit(1);
	}
	if(out_file=="") {
    int found=in_file.find_last_of('.');
		out_file=in_file.substr(0,found);
	}
	norm_vol * v;
	if(type=="uchar") {
    v=new vol<unsigned char>(dim_x,dim_y,dim_z);
	} else if(type=="char") {
    v=new vol<char>(dim_x,dim_y,dim_z);
	} else if(type=="uint") {
    v=new vol<unsigned int>(dim_x,dim_y,dim_z);
	} else if(type=="int") {
    v=new vol<int>(dim_x,dim_y,dim_z);
	} else if(type=="ushort") {
    v=new vol<unsigned short>(dim_x,dim_y,dim_z);
	} else if(type=="short") {
    v=new vol<short>(dim_x,dim_y,dim_z);
	} else if(type=="float") {
    v=new vol<float>(dim_x,dim_y,dim_z);
	} else if(type=="double") {
    v=new vol<double>(dim_x,dim_y,dim_z);
	} else {
	  cerr<<"Invalid type: "<<type<<endl;
	}
	v->load(in_file.c_str());
	norm_vol * trv=new tricubic_vol(*v);
	delete v;
	v=trv;
  lpt3d_tree tree;
	queue<lpt3d> q;
	for(int i=0;i<6;++i) q.push(lpt3d(i));
  while(!q.empty()) {
	  lpt3d s=q.front();
		q.pop();
		if(!tree.is_leaf(s)) continue;
		int level=s.orthant_level();
		if(level==max_level) continue;
		int m=min_level+6-level;
		if(m<3) m=3;
		if(level<min_level||(error(s,v,m)>max_error)) { 
		  tree.compat_bisect(s);
		  do {
		    q.push(tree.recent_code());
		  } while(tree.recent_next());
		}
	}
	tr1::unordered_map<point<3,double>,int> vtx_hash;
	vector<point<3,double> > vtx;
	int vtx_id=0;
	int tetra_count=0;
	tree.node_reset();
	do {
	  if(!tree.node_is_leaf()) continue;
		tetra_count++;
		array<point<3,double>,4> vs;
		tree.node_code().simplex(&vs[0][0]);
		for(int i=0;i<4;++i) {
		  if(vtx_hash.count(vs[i])==0) {
				vtx_hash[vs[i]]=vtx_id;
				vtx.push_back(vs[i]);
			  vtx_id++;
			}
		}
	} while(tree.node_next());
	ofstream prop_file;
	prop_file.open((out_file+".prop").c_str());
	for(int i=0;i<vtx.size();++i) {
	  prop_file<<((*v)(vtx[i][0],vtx[i][1],vtx[i][2]))<<" ";
	  prop_file<<(v->grad(vtx[i][0],vtx[i][1],vtx[i][2]))<<endl;
	}
	prop_file.close();
	ofstream node_file;
	node_file.open((out_file+".node").c_str());
	node_file<<vtx.size()<<" "<<3<<" "<<0<<" "<<0<<endl;
	for(int i=0;i<vtx.size();++i) {
	  node_file<<(i+1)<<" "<<vtx[i]<<endl;
	}
	node_file.close();
	ofstream ele_file;
	ele_file.open((out_file+".ele").c_str());
	ele_file<<tetra_count<<" "<<4<<" "<<0<<endl;
	int count=1;
	tree.node_reset();
	do {
	  if(!tree.node_is_leaf()) continue;
		array<point<3,double>,4> vs;
		tree.node_code().simplex(&vs[0][0]);
		if(tree.node_code().orientation()<0) swap(vs[2],vs[3]);
    ele_file<<count<<" ";
		for(int i=0;i<4;++i) {
			ele_file<<(vtx_hash[vs[i]]+1)<<" ";
		}
		ele_file<<endl;
		count++;
	} while(tree.node_next());
	ele_file.close();

	return 0;
}
